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(54) Method for correcting artefacts in a digital image 



(57) A method of correcting artefacts in a processed 
digital image represented by a multi-scale gradient rep- 
resentation being obtained by non-linear processing of 
a multi-scale gradient representation of an original im- 
age. A modified gradient representation is generated by 
modifying gradient images of the multi-scale gradient 



representation of the processed image so that the direc- 
tions of the gradient images are pixel-wise equal to 
those of the gradient images of the multi-scale gradient 
representation of the original image. A reconstruction al- 
gorithm is then applied to the modified gradient repre- 
sentation. 
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Description 

Field of the invention 

5 [0001] This invention relates to digital image processing. 

More particularly, it relates a method for correcting artefacts originating from non-linear treatment of components (co- 
efficients) of a multi-scale image representation. 

The method is applicable e.g. in a medical radiographic system, such as a computed radiography (CR) system or a 
computed tomography (CT) system. 

10 

Background of the invention 

[0002] In imaging systems where the final output image has to be evaluated by a human observer, a problem arises 
if the original image, as obtained from an image sensing device, contains detail information at various degrees of 

*5 coarseness, within a wide amplitude range. This situation is typical when the sensor has a good signal to noise ratio 
over a large dynamic range, which is the case with CR and CT In common working environments such as a hospital 
there is no time for selecting the optimal window/level setting (if any) for each individual image, so there is a need to 
display a single image on film or monitor screen, which reveals all the relevant diagnostic details with an acceptable 
contrast over the whole dynamic range. This problem has been recognized in the field of digital radiography, see: 

20 Maack I., Neitzel u., Optimized image processing for routine digital radiography, Proceedings International Symposium 
CAR "91. p. 109, Springer Veriag. 

[0003] Many methods for contrast enhancement have been developed, such as the commonly known techniques of 
unsharp masking (Cocklin M.L, Gourlay A.R., Jackson PH., Kaye G., Kerr I.H., and Lams P, Digital processing of 
chest radiographs, Imange and vision Computing, Vol. 1, No. 2, 67-78 (1983)) and adaptive histogram modification 
25 (Pizer S.M., Ambum E.P, Austin J.D., Cromartie R., Geselowitz A., Greer T, Ter Haar Romery B., Zimmerman J.B., 
and Zuiderveld K., Adaptive histogram equalisation and its variations, Computer vision, Graphics, and Image Process- 
ing, Vol. 39, 355-368 (1987)). These traditional methods suffer from two problems: 

1 . They fix to a particular filter size, whereas diagnostic details may occur at multiple scales in one and the same 
30 image. It may be desired to enhance details in the image at every possible scale. 

2. Artefacts are created in the vicinity of significant signal level transitions, e.g. at bone/soft tissue boundaries 
within the image. These artificats may impair the diagnosis. 

[0004] To overcome the first problem, multiscale methods for contrast enhancement have been developed in the 
35 past few years, see Vuylsteke P, Schoeters E., Multiscale Image Contrast Amplification (further referred to as MUSIC A 
algorithm), Proceedings of SPIE, vol. 2167, 5510560 (1994) and Lu J., Healy D M. Jr., Contrast enhancement via 
multiscale gradient transformation, In Proceedings of SPIE: Wavelet applications, Orlando, FL (1994); 
Fan J., Laine A., Contrast enhancement by multiscale and nonlinear operators, in VNfavelets in medicine and biology, 
Akram Aldroubi and Michael Unser eds., CRC Press (1996). 
40 [0005] While effectively curing this first problem, these methods also seem to perform well with regard to the second 
problem. The MUSICA algorithm is being used routinely in several hospitals around the world and there are no com- 
plaints about artefacts in CR. 

We have however noted that one has to be more careful in using these methods on CT images, since they have sharper 
signal level transitions and artefacts might still be created there. 

45 

Objects of the invention 

[0006] It is an object of the present invention to provide a method which corrects for artefacts originating from non- 
linear treatment of components (coefficients) in a multi-scale representation. 
50 [0007] Further objects of the present invention will become apparent from the description hereafter. 

Description of the invention 

Statement of the invention 

55 

[0008] The objects of the present invention are achieved by a method as set forth in claim 1. 

[0009] The method of this invention is based on a multiscale gradient representation such as described in: Mallat 
S., Zhong S., 
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Characterization of signals from multiscale edges, IEEE Trans, on Pattern analysis and Machine Intelligence, Vol. 14, 
No. 7 (1992)). 

[0010] The method has two basic steps. The first step is the nonlinear processing step, e.g. a contrast enhancement 
step (e.g. using an algorithm similar to that of Lu and Healy (see above)). 

The second and novel part is an artefact correcting step. This step can be used to correct artefacts caused by any 
nonlinear treatment of the coefficients in the multiscale gradient representation and its use is therefore not necessarily 
restricted to the problem of contrast enhancement. 

Detailed description of the present invention 

Multiscale gradient representation 



[0011] A digital image f is a map from a square grid of pixels Pto the real numbers R . The multiscale gradient rep- 
resentation is a decomposition of the image into a series (say K, K=1 ,2,3, ...) of vector-valued images on P(i.e. maps 
is from PtoR 2 ) W 0 f, and a residual image S K f , each having the size of the original image (no subsampling 
takes place). Let SqI- /and for j= 1, ... , K, Sf \s obtained from S^f by applying a smoothing filter, then W/is defined 
to be the (discrete) gradient of S/ . We will refer to W 0 f , .... V& K . A /, S K f as the multiscale gradient representation (MGR) 
of digital image 

[0012] The smoothing and gradient filtering are obtained by two-dimensional convolution. 
20 An efficient implementation of such filters is based on successively applying a horizontal 1 -D convolution filter, followed 
by a vertical 1-D convolution filter, applied to the result of the first filter. If the horizontal filter is denoted by a and the 
vertical filter by b, then the resulting 2-D filter is commonly denoted by a®6. 

[001 3] The filters at multiple scales are generated from a basic filter at the smallest scale. The corresponding filters 
at subsequent larger scales j=1, ...K, are obtained by inserting 2>-1 zero-valued coefficients between any two coeffi- 
zs cients of the basic filter. 

[0014] In a preferred embodiment, a basic smoothing filter h 0 is used with coefficients (1/16, 4/16, 6/16 , 4/16, 1/16). 
The underlined filter tap is defined to be the filter origin (i.e. center tap, which has index zero). The basic gradient filter 
is chosen as g n =f1/2,-1/2). 

The corresponding filters at subsequent larger scales are obtained by inserting zero coefficients as described above: 

30 

h, = (1/16,0,4/16,0,6/16,0,4/16,0,1/16), = (1/2,0,-1/2), 
h 2 = (1/16,0,0,0,4/16,0,0,0,6/16,0,0,0,4/16,0,0.0,1/16). 
g 2 = (1/2,0,0,0,-1/2). 

35 if a filter has an even number of taps (e.g. c/ 0 ), then the filter at the next larger scale has an odd number of taps, and 
the filter origin is shifted one position to the right. 

Appropriate combinations of these filters in horizontal and vertical directions are selected to generate the 2-D filters 
needed for the MGR decomposition: 



S 0 f=f 



S/= (h^ ®n Al ) * S. 1 f , for 0<j<K 
W/= «g y ®8 0 rs/, (5 0 ®g y ys/) (1) 



Here, * is the convolution operator and 5 0 the Kronecker delta (i.e. the 1 -D filter in which the center tap has a coefficient 
equal to 1 and all other taps have a coefficient equal to 0). 

Hence, g y <8>5 0 means filtering by p y in the horizontal direction only and the 6 0 ®gy is equivalent to applying the filter g } in 
so the vertical direction only. 

[0015] Clearly S/is a smoothed version of S^f , while T^/is the (discrete) gradient of S/aX scale j. 
When convolving an image f with a filter, the image is extended at the borders in a continuous way by reflection about 
the image borders. After performing the convolution with this extended image, the result is again restricted to the 
pixelset of the original image. 

ss [0016] The decomposition is schematized in Figure 1 ( in this Figure, superscripts (1 ) and (2) refer to the horizontal 
and vertical component of the gradients respectively). 

It is well known that this decomposition is reversible, i.e. the original image can be entirely recovered from the repre- 
sentation by a reconstruction procedure, denoted by the inverse multiscale gradient (IMG). To this end, two additional 
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filters p 0 and q 0 are needed with the constraint that 

P 0 * + % * 9 0 = 5 0 (2) 

5 

p Q and g 0 are basic 1-D filters. 

[0017] Filters p ; and q } at larger scales j are generated by inserting 2M zero coefficients between successive taps, 
as described above. 

[0018] The IMG reconstruction is carried out as explained in Zong, S, "Edge representation from wavelet transform 
10 maxima", Ph.D. Thesis, New York University, 1990 as follows: 

S/ = (p. @ p.) * S y+1 f + (q. 9 /■) * IV/ 1 V + (/ y 0 q.) * w} 2) f, j>0. (3) 

15 The 1-D filter l y is defined as ^* - 1/2(8 0 + Pfty. (superscripts (1) and (2) denote horizontal and vertical direction). 

[0019] The image is recovered from the low-resolution image S^f and the gradients I# 0 /,...,W K _ 1 /by repeated use 
of this equation for j=K-1 ,...,0. (see Figure 2). 

[0020] The filters p 0 and q 0 should be sufficiently smooth to avoid discontinuities when reconstructing the images 
from a set of nonlinearly modified transform coefficients (as will always be the case when doing contrast enhancement). 
[0021] In a preferred embodiment p 0 = (1/4,0,1/2,0/1/4) and 1/32, -5/32,- 13/32, -25/32, 25/32 , 13/32,5/32, 1/32). 
These filters comply with constraint (2). 

[0022] I n a first embodiment, the number of pixels remains constant at each scale. Such a representation is commonly 
called a multiscale stack. 

[0023] In a second embodiment, the numbers of pixels is lowered in each step of the decomposition, by subsampling. 
25 Such a representation is commonly denoted by a muttiscale pyramid. 
[0024] The MGR decomposition is then defined as : 
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S 0 f = f 

Sjf = l(h 0 ®h 0 )*S j _J y for l<o<K (4> 

W// = ((£ 0 ®5 0 )*S ; /,(S 0 ®g 0 )*£./) for OSlSjSk-l 



[0025] Here I denotes downsampling by a factor of two, keeping only the pixels with even coordinates (2x,2y), where 
(x,y) refers to the pixel coordinates in the current grid. 

[0026] The pixels of the resulting image are renamed from (2x,2y) to (x,y). The number of pixels is halved in both 
dimensions, in each stage of the pyramidal decomposition. The reconstruction from the MGR pyramid is done using 
40 the equation 

iS/= (p 0 ® p 0 )*S M + Mq 0 * / o r^ (1 V + l(! Q <8> q 0 yw} 2) f, j>>0 (5) 



where Pq = (1/4,1/2,1/4). 

[0027] The reconstruction procedure (5) is equivalent to subsampling the result of the reconstruction (3), taking into 
account that p 0 <8> p Q ) -l 1- (p 0 <s> p Q )f, for any arbitrary image I if p 0 has zero coefficients between subsequent taps. 
The latter condition is fulfilled if Pq = (1/4,0,1,2,0,1/4) 

[0028] Equation (5) tells us how to obtain a subsampled version iSyfof SyA. The full-size reconstructed image £>.f is 
next calculated from island Wj{Wjfbe'\nQ the gradient of SyO in the following way. 

First the pixels with coordinates (x,y) of the subsampled reconstruction is\f (further called sample pixels) are assigned 
new coordinates values (2x,2y). 

[0029] The pixels of the full-size reconstructed image at positions (2x+1 ,2y) lying in between the two sample pixels 
(2x,2y) arid (2x+2,2y) are computed by: 
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S]f(2x + l.2y> = }/2[s)f(2x,2y) -2W j m f(2x,2y) + S]f(2x + 2,2 y) + 2 f{2x + \,2y)] ( 6 ) 



5 [0030] The pixels with coordinates (2x,2y+1 ) lying in between the two sample pixels (2x,2y) and (2x,2y+2) are com- 
puted by 

S)f(2x.2y+ \)=l/2[Sjf(2x2y)-2W i ,2 >f(2x2y) + S)f{2x,2y + 2) + ^^ (7 ) 

[0031] The pixels with coordinates (2x + 1 ,2y + 1 ) in between the four sample pbcels(2x,2y),(2x+2 l 2y),(2x,2y+2) and 
(2x+2,2y*2) arc computed by: 
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5*/(2x + l.2y^l) = l/4 



S)f(2x + l,2y) - 2Wj (2) f(2x + 1,2 y) + 
S,/(2x,2y + 1) - 2^: (,> /(2x,2y + I) + 
S 7 /(2x 4- 1,2 y + 2) + 2 W / : <2) /(2.t + L2y + 1) + 
S]f(2x + 2,2y + 1) + 2W: <,) /(2a* + L2.y + 1) 



(8) 



25 [0032] Wo will rclcr to the latter procedure (described by equations 6-8) as gradient based interpolation (GBI). 
[0033] The pyramidal reconstruction scheme is depicted in Figure 4. 

Contrast enhancement step 

30 [0034] First the mulliscale gradient representation "^ 0 f,...,'^ K . 1 f,S K /(or W Q f, ... t V&^ f,S' K f\n case of the pyramid) of 
the image / is calculated. 

[0035] To enhance the contrast in the image we want to modify the transform coefficients in such a way that details 
of small amplitude are enlarged at the expense of the ones with larger amplitude, and this uniformly at all scales. 
Therefore the transform coefficients must be modified in a nonlinear, monotonic way. 
35 [0036] In preferred embodiment, the following contrast enhancement function is used: 



40 



,-lmY 



if 0<x<c 



6 J 



(9) 



m is an upper bound for the transform coefficients,0<p<1 determines the degree of nonlinearity and c prevents the 
45 function yfrom getting too large for small values of x. This is done to prevent unnecessary over-enhancement of the 
noise in the image. The value of c should therefore be related to the noise content of the image, e.g. we can take cto 
be the estimated standard deviation of the noise in the image. 

[0037] The function y is continuous and monotonically decreasing, has a constant value 



so 



for 0<x<cand reaches the value 1 for x=m. (see Figure 5). 

[0038] The modified multiscale representation "v^,..., "P f K . 1 , S K f (alternatively "5^,.. ^ K _^ ,S' K f \n case of a pyramidal 
decomposition) is defined as 
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v/0 = ^/(0>([(wf /(<)) 2 + (w^foifj 2 ) do) 

5 Vj (o = Wjf(o.y([(w; (,) f(of +(^: (2) /(i)) 2 ]" 2 ) 

for all pixels / and all scales 0<j<k 

[0039] This has the following effect: gradient vectors with small magnitude are multiplied by a larger factor than 
10 gradient vectors with large magnitude, so details of small amplitude are favoured at the expense of more striking ones. 
[0040] The enhanced output image f is then calculated by carrying out the IMG on "v^, V KV Srf (or , 

If necessary, the grey values are then linearly rescaled to the original range. 

[0041] Note that in this modification, only the gradient magnitudes are affected, not their directions. This is done to 
15 prevent artefacts upon reconstruction. However, this goal is only partly achieved, because the method treats the trans- 
form coefficients as if they were independent of each other. 

However this is not entirely true, since the gradient at the y-th scale is in fact a smoothed (possibly subsampled) version 
of the gradient at the /-1-th scale. If the gradients are independently modified in the above way, the latter property is 
lost, even more, V r 0t ... t 'P K ^ t Sf(f\s no longer a MGR. 
20 This means that there exists no image g such that V 0 ,. „ty K _^S K f is the MGR of g. 
The same holds of course in the case of pyramidal decomposition ^ 0 ,- -,^ .S^f*. 

As a consequence, the gradients of the reconstructed image /will no longer be parallel to those of the original image 
fas we envisaged. 

[0042] The distortion of edge gradient orientation due to independent modification of MGR pixels at subsequent 
25 scales, is an unwanted effect, and may give rise to visible artefacts. 

Artefact correction step 

[0043] The above problem is solved if the reconstruction is applied to a set of gradient images ~% Q !% A ,...^ k .^ instead 
30 of ~v* Qt . . . J&k-i (or V Q .... , ) which comply to the following constraints 

D1 . there exists an image gsuch that ^ 0 ,^ 1 ,. .,^ /M( S K g is the MGR of g 

02. 



35 



40 



XjiO^WjfO) IHfrj/iO) for a11 0<j<K and all pixels i 

D3. ~%j is more equalized than V&f{Wjf). 



[0044] In D3, the term 'more equalized 1 is to be understood in the sense of equation (10): ~k } is called more equalized 
than PPyfwhen, if 

is 'small' then Lfy/)l>lT^/)l, while 
if lW/(/)l is 'large', then l^//)l<lW/(/)l. 
45 [0045] In other words, the artefact correction method aims at redirecting (D2) the gradient vectors of the enhanced 
(D3) image to the gradient directions of the original image in a way consistent with the multiscale gradient representation 
structure (D1). 

[0046] We describe below two alternative embodiments of methods that construct such a set ^ 0 »^1 •■ ->%k-v 

50 - Alternating convex projections, 

[0047] It is well known that projecting alternatively onto two convex nondisjoint subsets of a linear space is a con- 
vergent scheme that yields a point in the intersection of these two sets (see Youla D., Webb H., Image restoration by 
the method of convex projections, IEEE Trans, on Medical Imaging 1, 81-101 (1982)). 
55 [0048] This knowledge is used in a first embodiment to construct gradient images satisfying the above demands 
D1.D2 and D3. 

[0049] The method starts from the modified representation, V Ql .. ^ K .^ t S K f {ox from "0^,. ..."v^ .S> K f) defined in (10) 
and iterates an alternating sequence of projection operations onto the two vector spaces as follows : 
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1. projection onto the space W of multiscale gradients, i.e. all sets of vectorimages 

such that there exists an image g for which ^ y = T&jg (or ~% s = ^j9) for 0<j<K 

2. projection onto the space V|| of sets of vector images 

w 




is such that for all pixels i and all scales 

0<j<K t X^Wjfii) or ( ^.(i)|^;/(0 ) . 

20 

[0050] The result of the first projection is fed into the second projection. The resulting vector-images are again fed 
into the first projection stage, and so on. 

[0051] The first projection onto the space W of multiscale gradients is done by performing on a set of vector-images 
the multiscale gradient reconstruction procedure (equations (3) or (4), (5), (6), (7) and subsequently applying to the 
2S resulting image the multiscale gradient decomposition procedure (equation (1) or (4)). This projection operation is 
denoted by P w 

[0052] The second projection operation onto the vector space of images of which the vector-valued pixel values 
are parallel to the corresponding gradient pixels of the multiscale gradient representation Wf t is done by taking the 
orthogonal projection onto Tftfli) (Wjf{i) for every /and / The orthogonal projection of a vector v onto the vector w is 
30 equal to 
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Projection onto V| t , is denoted by P^. 

[0053] The output image is obtained by applying the IMG reconstruction process to the result of iterating the above 
alternation of projection operations 

P vl °P w ..P vl 0 P w {V Q ,...,V K _„S K f} (11) 

and analogous in case of the pyramid scheme. 
45 [0054] We need approximately 10 iterations for convergence, i.e. we need to do the operations P w and P^ 10 times. 
This algorithm results in a multiscale gradient representation that satisfies all three above constraints D1 .,02. and D3. 

- Multiscale gradient projection 

so [0055] A second embodiment of the artefact correction process is described next, which is based on the multiscale 
stack embodiment of the MGR; subsequently a third embodiment of artefact correction based on the pyramidal repre- 
sentation will be disclosed. 

[0056] The second embodiment is based on the fact that for every /, 0 </< /C-1 there exists a linear operator f^on 
the space of vectorimages such that for every image f 
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This means that the multiscale gradient image at a specific scale can be computed from the multiscale gradient image 
at the preceding smaller scale by a linear filter operation Fj[) 

Indeed, 

[0057] 

%xf = (<£; +I ®So)*S,. +l /. (8 0 ®g >+1 )*S J . l f) 

= <(£; + ,®5o)*(A,®fc;>* V- tfo®g J+l )*(h j ®h J )*S J f) 

= ( ( *y, * hj&h, ) * Sjf , ( hjVg^ *h,)* SJ ) 

= ((b j *g J ®h j )*S j f , (h j <8>b j + g j )*S J f 
= Ub j ®h j )*lg j ®S 0 )*S j f , (h j <8>b j )*{8 0 ®g j )*S j f ) 
= {(b^h,), (h j ®b j ))*Wjf (13) 

Fj thus acts on a vector image ^ as 

Ffft = (bj*h Jt hj*bj)*V (14) 

The result of this operation is also a vector image. 
From (13), the filter bj is defined by the relation 

Y9j=9u,* h i < 15 > 

Hence the filter coefficients of 6 y can be calculated from h$ and g 0 . This is done in the following way. 
55 Filters are often represented by their z-transform. The z-transform A(z) of a filter a is a function of the complex number 
z defined by A(z) = Za k z k , where the sum runs over all fitter taps k, in which a k are the coefficients of the filter kernel. 
With our choices of ho and g 0 , 

40 H Q {z) = in6(z~ 2 +4z^-&+4z+z 2 ) and G 0 (z) = 1/2(1-z) 

[0058] One of the advantages of the z-transform is that it converts a convolution into a product, i.e. for two filters a 
and b the z-transform of a*b is A(z)B(z). 

Therefore bj can be easily calculated by means of its z-transform. From (15) 
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B y (z)=G y+1 (z)H,(zyG y (z). 



[0059] This yields 



B Q (z)=1/16 (z- 3 +5z- 2 +10r 1 +10+5z+z 2 ) t thus 
b 0 = (1/16,5/16,10/16,10/16,5/16,1/16), 

= (1/16,0,5/15,0,10/16,0,10/16,0,5/16,0,1/16), etc.. 

[0060] The artefact correction method according to a second embodiment starts by decomposing the enhanced 
image f into its multiscale gradient representation, yielding l5 0 ,..., ^?k-i»^k* 

[0061] In accordance with the criterion D2, these vectors have to be redirected in the direction of the multiscale 
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gradients of the original image. 

[0062] For /=0,. .../C-1 , let Pj denote the pixelwise orthogonal projection operator onto V3j. 1 -Pj is then the pixelwise 
orthogonal rejection w.r.t. V&f. The orthogonal rejection v ± of a vector v w.r.t. the vector w is given by the formula 

5 . . v°V 1) + v (2 V 2) - 

[0063] As a resutt of projecting 75 0 onto W 0 f according to P 0 some information is lost and this can lead to a loss of 
10 contrast upon reconstruction. 

This lost information resides in (1-Po)A) 



if us rui>i iinuiiiiciuun itr^iuu^ in (itq;uq. 

(1 - P 0 )V 0 has zero component along W Q f. 
However, 



is 



does have a nonzero component along W 0 f. 
And, 

20 



although it belongs formally to scale 1 , can be considered as a contribution to scale 0. since it is in no way a smoothed 
2S version of P 0 U 0 (it is obtained from 0-P 0 )Tf 0 ). 
Therefore we add 



(\>2)P 0 F 0 {(\-P 0 )0 0 } to P 0 U 0 



30 



35 



40 



[0064] The factor 1/2 comes in because the filter b 0 used in the operation F 0 has gain (=sum of all filter taps) equal 
to 2, so by adding just 



we would add a contribution that is 2 times too large. 

[0065] By repeating this procedure the vector image at smallest scale ^ 0 is defined by: 



2 4 

45 1 

2 



+ 7^PoF K -20-Po)...F 0 {\-P 0 )U 0 



(16) 



[0066] This will be the full contribution to scale 0. 
so The vector image at the next scale "Xy is obtained in a simular way, starting from 



Fo{*o} 
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(17) 



10 



15 



and so on, and similarly 

**2 = PK^ F K.s *K-3) + J P*_ 2 F K _ 2 (1 - P K . 2 )(F K _ 3 ^ 3 ) (18) 
^-1 = ^ ( F K. 2 ^K. 2 ) 09) 

[0067] The final output image is then reconstructed by applying the IMG to 

[0068] According to our findings, 

2$ {x 0 , . . . , X K _i , S K J- 

satisfies constraints D2. and D3. exactly and D1. approximately 

[0069] This, algorithm is faster than the alternating convex projection algorithm and preserves better the contrast 
gained in the contrast enhancement step. 
30 [0070] We now indicate what changes should be made when using the pyramid decomposition scheme, according 
to a third embodiment of artefact correction. 

[0071] There still exists a linear operator £ Q such that for every image f 

3S %i'=F o 0^0 (20) 

but in this case /=* 0 acts on a vectorimage ~P as 

40 F 0 {V)= l(b 0 ®h 0 ,h 0 <8>b Q )*~& (21) 

[0072] The filter b Q is the same as defined before, but the operator P Q also involves subsampling, denoted by i. The 
discarded components of projection (1-Py) are added back in a way similar to the second embodiment, but the effect 
of subsampling has to be encompassed by proper interpolation. 
45 [0073] In our preferred embodiment bilinear interpolation was chosen, denoted by the operator /. This operator is 
applied to each of the contributions at larger scales, in order the match the image dimension of the next smaller scale. 
[0074] The corrected gradient images Xj are then obtained by : 

so % = Po&o + \ Po^oO - P 0 )tk + \ ^(AoO - ^o) F o( 1 * 

+ ^i P o /Ar1/r oO " Po)-W " p o)A>) 
55 % - Pi(F 0 %) + l P,IF 0 V - P,)(F 0 X) + \ p/f 0 V - ^)F 0 (1 - P,)(F 0 %)+... 
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**-2 = F K . 2 (F 0 ^. 3 ) + |P^ 2 /F 0 (1 - Pk_ 2 )F 0 * k _ 3 ) 
= P/c^Fo^/c-a) 



(22) 



[0075] Of course, when projecting a gradient field at some larger scale onto a gradient field at a smaller scale, we 
have to identify the corresponding pixels, i.e. if ¥ is Ac scales larger than 2" the pixel (x,y) of T' should be identified with 
the pixel (x2 k ,y2 k ) of 1 . 



1. A method of correcting artefacts in a processed digital image / represented by a multi-scale gradient representation 



comprising gradient images at multiple resolution levels and a residual image, said gradient representation being 
obtained by non-linear processing of a multi-scale gradient representation 



of an original image f , 
characterised by the steps of 

(i) generating a modified gradient representation by modifying the gradient images of Vto make their directions 
pixel wise equal to those of W, and 

(ii) reconstructing a final output image from this modified gradient representation by applying an inverse mul- 
tiscale gradient transform (IMG), said inverse transform being such that it yields / when being applied to W. 

2. A method according to claim 1 wherein 

a multiscale gradient representation 



is generated and wherein said final output image is reconstructed from said multiscale gradient representation 



Claims 



w={w 0 /,...,w ir . 1 /,s Jf / } 




said multiscale gradient representation 




satisfying the following constraints: 



(1 ) an image g * f exists such that for/= 0,...,K- 1 , represents the gradient image of pat the y-th scale, 

(2) in every pixel ^-has the same direction as the multiscale gradient image Wjfol the original image fat 
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the /-th scale, 

(3) if the original multiscale gradient image Wj has small magnitude, then the magnitude of ^ y is larger 
than the magnitude of the original gradient image; 

if the original gradient image has large magnitude, then the magnitude of 1( } is smaller than the magnitude of the 
original gradient image, at each pixel and at each scale j. 

3. A method according to claims 1 or 2, wherein the multiscale gradient representation of the original image is obtained 
by (i)applying one-dimensional gradient filters independently to the rows and columns of a digital image represen- 
tation of the original image thereby yielding the horizontal and vertical components of the muttiscale gradient image 
at the smallest scale, and applying a two-dimensional low-pass filter to said digital representation of the original 
image,thereby yielding an approximation of the original image at a next larger scale, (ii)identifying the above op- 
erations as the first step of an iterative loop, and performing additional iterations using the approximation image 
resulting from a previous iteration instead of the original image, with the filters being adapted to the current scale, 
thereby yielding gradient and approximation images of the original image at the successively larger scales. 

4. A method according to claim 3, where said gradient filters have filter coefficients (0.5 , -0.5) or a multiple thereof 
and where the low-pass filtering is obtained by applying a horizontal one-dimensional low-pass filter to the rows 
of an image, followed by a vertical one-dimensional low-pass filter applied to the column ns of the resulting image 
of the horizontal filter, both low-pass filters having coefficients selected from the following set: (0.5 , 0.5), (0.25, 0.5, 
0.25), (0.0625, 0.25, 0.375, 0.25, 0.0625). 

5. A method according to claim 2, wherein constraints (1 ) to (3) are fulfilled by 

(A) obtaining an initial estimate of the multiscale gradient representation X by: 

a) applying an inverse multiscale gradient transform (IMG) to the multiscale representation V of the proc- 
essed image /, and subsequently decomposing the resulting image into a muttiscale gradient represen- 
tation X, 

b) redirecting all vector-valued pixels of X obtained in (a) by pixelwise orthogonal projection onto the pixels 
of the multiscale gradient representation W of the original image, yielding a multiscale grandient repre- 
sentation X", so that the pixels of X and W have identical directions at all scales, 

(B) identifying the above operations as the first step of an iterative loop and performing additional iterations 
using the result A" of step (b) from the previous iteration as the input of step (a) of the current iteration, instead 
of V, thereby yielding subsequent improved estimates of the multiscale gradient representation X , 

(C) obtaining a multiscale gradient representation Xas the result X* of step (2), after the last iteration. 

6. A method according to claim 2 wherein constraints (1 ) to (3) are fulfilled by : 

applying an inverse multiscale gradient transform (IMG) to the multiscale representation Vof the processed 
image / and subsequently decomposing the resulting image into a multiscale gradient representation 



obtaining ^ 0 as the sum 




where T 0jt is obtained from t) 0 by applying the following steps: 
Take T) 0 to be the first input of A1 and 
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A. Repeating for j=0 to k: 

1. extract the component of the input vector image parallel to W 0 f and the component of the input 
vector image orthogonal to T^ 0 f by pixelwise orthogonal projection of this input vector image onto 
ffl Q f, the parallel component being the output vector image of A1, 

2. if j<k, apply an operator Fj to the orthogonal component, the operator Fj being such that it yields 
Wj+i f when acting on Wf, divide this result by a factor of 2, use this result as the next input to step A1 , 

B. Taking the last output of A1 , corresponding with loop index j=k, and if necessary interpolating this result 
to get the same number of pixels as Wrf, the output of B being T 0k , 

for j=1,.„, K-1, obtaining ^ as the sum 



is - - 

where is obtained from F M ) by the following procedure : 
20 Take (F M (X M ) to be the first input of C1 

C. Repeating for j=i to K: 

1. extract the component of the input vector image parallel to *?/and the component of the input 
25 vector image orthogonal to T^/by pixelwise orthogonal projection of this input vector image onto W-f t 

the parallel component being the output to C1 ( 

2. If j<k, apply an operator Fj to the orthogonal component, 

the operator F ; being such that it yields "ffi^f when acting on W/and divide this result by 2, use this 
30 resulting vector image as the next input of C1 . 

D. Taking the last output of C1 , and if necessary interpolate this result to get the same number of pixels 
as V&f. 

7. A method according to claim 1 wherein said non-linear processing is a contrast enhancing processing. 

35 
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